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Abstract 

We perform numerical optimization of the axisymmetric flows in a sphere to minimize the critical 
magnetic Reynolds number Rm cr required for dynamo onset. The optimization is done for the 
class of laminar incompressible flows of von Karman type satisfying the steady-state Navier-Stokes 
equation. Such flows are determined by equatorially antisymmetric profiles of driving azimuthal 
(toroidal) velocity specified at the spherical boundary. The model is relevant to the Madison 
plasma dynamo experiment (MPDX), whose spherical boundary is capable of differential driving 
of plasma in the azimuthal direction. We show that the dynamo onset in this system depends 
strongly on details of the driving velocity profile and the fluid Reynolds number Re. It is found 
that the overall lowest Rm cr « 200 is achieved at Re « 240 for the flow, which is hydrodynamically 
marginally stable. We also show that the optimized flows can sustain dynamos only in the range 
Rm cr < Rra < Rm C r2, where Rm cr 2 is the second critical magnetic Reynolds number, above which 
the dynamo is quenched. Samples of the optimized flows and the corresponding dynamo fields are 
presented. 



I. INTRODUCTION 



The creation of specific flows of an electrically conducting fluid is the key element and ma- 
jor challenge in experimental investigation of dynamo phenomenon. Appropriate candidates 
for this role are the flows capable of the dynamo action, a number of them are theoretically 
studied in literature [1—10] . In most of these studies kinematic dynamos are considered, 
i.e., the magnetic induction equation is treated as an eigenvalue problem for an unknown 
magnetic field and a prescribed laminar flow. The model flow is usually chosen in a simple 
analytical form, which is not necessarily determined from the fluid dynamics. In fact, the 
majority of analyzed flows leading to kinematic dynamos do not satisfy the Navier-Stokes 
equation (e.g., flows from Refs. [1, 3-6, 9]), so their structure cannot be reproduced exactly 
in the laboratory. Nevertheless, it is possible to obtain experimentally a flow sufficiently 
close to the model one, which may result in the dynamo action. 

In the past decade, several groups constructed dynamo experiments intended to achieve 
such flows with liquid metals [11-17]. Although the obtained flows were highly turbulent 
in all the experiments, their mean parts were expected to sustain a dynamo field. However, 
only three experiments were successful in dynamo demonstration: experiments in Riga, 
Latvia [11], and Karlsruhe, Germany [13], where the flows were strongly constrained and 
the influence of turbulence was small, and a von Karman sodium experiment in Cadarache, 
France, where ferromagnetic impellers played the critical role [16, 18]. It seems likely that 
hydrodynamic turbulence in unconstrained flows significantly inhibits dynamo onset. 

The presence of turbulence is an inevitable problem in all liquid metal dynamo exper- 
iments. This is due to the extremely low magnetic Prandtl numbers of the liquid metals, 
i.e., the ratio of kinetic viscosity v to resistivity i] or, equivalently, the ratio of magnetic 
Reynolds number to fluid Reynolds number Pm = vjr\ = Rm/Re (e.g., for liquid sodium 
Pm ~ 10~ 5 ). In order to reach the magnetic Reynolds numbers sufficient for dynamo excita- 
tion (Rm ~ 10 1 — 10 2 ), very high fluid Reynolds numbers (Re ~ 10 6 — 10 7 ) are required. As 
a result, the corresponding flows in experiments are always turbulent, making it difficult to 
achieve dynamos and to compare experimental data with predictions of laminar kinematic 
theory. Note that the liquid metal laboratory dynamos with Pm <C 1 are in the regime of 
the solar convection zone and the interiors of planets, while hot accretion disks and galaxies 
have Pm ^> 1. 
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The present study is motivated by construction of the Madison plasma dynamo experi- 
ment (MPDX, Fig. 1), which is designed to investigate dynamos excited by controllable flows 
of plasma [10, 19]. The use of plasma as the electrically conducting fluid gives experimental- 
ists flexibility in choosing a regime of dynamo operation. By adjusting experimental controls, 
one can change driving velocity, density, electron and ion temperatures, etc. This makes it 
possible to adjust Pm, Rm and Re at will and study laminar dynamos with Pm ~ 1 and 
Rm ~ Re ~ 10 2 (Table I). Such flexibility is advantageous over the liquid metal dynamo 
experiments. The experimental vessel in MPDX is an aluminum sphere of 3 meter in di- 
ameter [Fig. 1(a)]. Plasma is confined by an axisymmetric multicusp magnetic field created 
by 36 equally spaced rings of permanent magnets with alternating polarity. The plasma 
filling the vessel is mostly unmagnetized since the multicusp field is localized near the vessel 
wall. The novel feature of the experiment is the mechanism for creating controllable plasma 
flows [Fig. 1(b)]. An electric field applied across the multicusp magnetic field drives the 
edge of the plasma azimuthally, while viscosity couples momentum from the edge to the 
unmagnetized core. Nearly arbitrary profiles of plasma azimuthal (toroidal) velocity v^(9) 
can be imposed at the sphere's boundary by modulating the electric field as a function of 
polar angle using discrete electrodes. This concept of plasma stirring has been successfully 
tested in the plasma Couette experiment (PCX) [21], and it allows a unique way to conduct 
laboratory studies of various astrophysical phenomena including the dynamo [10, 22], the 
magnetorotational instability [23], and the Parker instability [24]. 

As shown in Refs. [25, 26], the toroidal motion alone does not sustain a dynamo magnetic 
field, some poloidal flow is necessary. In a bounded sphere the poloidal flow develops self- 
consistently from differential toroidal rotation, which is controlled in MPDX by boundary 
driving velocity v$(Q). It can be shown using the Navier-Stokes equation that the intensity 
of poloidal flow is determined by fluid Reynolds number Re: the poloidal flow is normally 
stronger for the larger values of Re. There is a minimum amount of poloidal motion necessary 
for dynamo action [27]. From this point of view, the large values of fluid Reynolds number 
Re are desirable in the experiment. At the same time, to avoid turbulent regime of the flow, 
the values of Re should be less than the hydrodynamic instability threshold. We emphasize 
here again that fluid Reynolds number Re in the unmagnetized MPDX plasma can be varied 
by changing ion parameters of plasma: density n , temperature Tj and mass ^. 

The poloidal motion is required for dynamo excitation, but does not guarantee it. The 
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FIG. 1: The Madison plasma dynamo experiment (MPDX): (a) a sketch of the experiment; (b) 
the electrode configuration near the wall for driving plasma velocity v^(6). The spherical system 
of coordinates (r, 9, 4>) is shown. Center line (CL) corresponds to the axis of symmetry. Reprinted 
with permission from Phys. Plasmas 19, 104501 (2012) [33]. Copyright 2012 American Institute 
of Physics. 

dynamo onset is very sensitive to the details of flow structure and corresponding profile 
of driving azimuthal velocity v^d). The ability to create arbitrary profiles of v^O) in 
MPDX raises the question of their optimization for the dynamo excitation. The goal of 
our study is to find numerically the optimized profiles of driving velocity v${9) and the 
corresponding equilibrium flow structures in the model relevant to MPDX. We refer to a 
flow as optimized if it minimizes the critical magnetic Reynolds number Rm cr required for 
the kinematic dynamo onset. Physically, a lower Rm cr means a lower driving velocity and/or 
lower electron temperature (higher resistivity) of the plasma needed for achieving dynamo. 

The problem of flow optimization for the kinematic dynamo in spherical geometry was 
addressed by many researchers [28-32]. The necessity of balancing the relative amplitudes 
of the toroidal and poloidal flow components for the dynamo action was originally noticed 
for the Kumar and Roberts flow (KR flow) [3] and later for the Dudley and James flows 
(DJ flows) [4]. More rigorously, the idea of flow optimization was introduced by Love and 
Gubbins [28, 29], who optimized the relative amplitudes of the four spherical harmonic 
components of the KR flow, keeping their radial dependences fixed. O'Connell with co- 
authors [30] and Holme [31] optimized DJ flows, allowing the radial structure of spherical 
harmonic components to vary. Gubbins and co-authors showed that the optimized flow of 
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TABLE I: Expected parameters of MPDX. Dimensionless numbers Re, Rm and Pm are estimated 
from the Braginskii equations [20] (see corresponding formulas in Refs. [10, 22]). 



Quantity 


Symbol 


Value 


Unit 


Radius of sphere 


Ro 


1.5 


m 


Peak driving velocity 


v 


0-20 


km/s 


Average number density 


n 


10" _ 10*9 


m ° 


Electron temperature 


T e 


2-10 


eV 


Ion temperature 


Ti 


0.5-4 


eV 


Ion species 




H, He, Ne, Ar 




Ion mass 


Hi 


1, 4, 20, 40 


amu 


Fluid Reynolds 


Re 


- 10 5 




Magnetic Reynolds 


Rm 


- 2 x 10 3 




Magnetic Prandtl 


Pm 


10~ 3 - 5 x 10 3 





KR type maximizes the non-axisymmetric part of kinetic helicity (see Ref. [32] and references 
therein). 

It appears that all flow optimizations considered in the dynamo literature so far are per- 
formed for simplified model flows of a particular type (either KR or DJ flows). The structure 
of these flows is usually prescribed by their type and is not determined self-consistently from 
the hydrodynamic equations. In this sense, the optimized flows found are not realistic and 
cannot be reproduced in an experiment. In contrast to these studies, in the present paper we 
deal only with realistic flows found as solutions to the steady-state Navier-Stokes equation 
with boundary conditions specified by the driving velocity profiles v^{6). 

The structure of the paper is as follows. In Sec. II, we briefly describe the model of MPDX 
used in our study (the basic equations of the model and methods of their solution are given 
in Appendix A). In Sec. Ill, we discuss the factors that influence the dynamo threshold. In 
Sec. IV, the results of the flow optimization are reported. In Sec. V, we summarize our main 
findings. 
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II. MODEL AND OPTIMIZATION PROCEDURE 



We perform our study in the framework of single-fluid magnetohydrodynamics (MHD), 
which is a good approximation for the MPDX plasma. For simplicity, we ignore the effects 
of plasma compressibility and the details of plasma confinement and drive near the wall. 
We neglect the multicusp magnetic field and the applied electric field and assume that the 
velocity profile is specified at the sphere's boundary. As shown in Ref. [22] for the model 
relevant to PCX (a cylindrical prototype for MPDX), these ignored details only play a role 
in the relatively thin boundary layers, and are not essential in the bulk of unmagnetized 
plasma. 

The equations of our model in non-dimensional form are 

= V 2 v-i2e[(v- V)v + Vp], (1) 

7„v = V 2 v - Re [(v • V) v + (v • V) v + Vp] , (2) 

7bB = V 2 B + RmV x (v x B), (3) 

= V • v = V • v = V • B, (4) 

where v = V/Vo and p = P/(p Vq) are normalized velocity and plasma pressure in equi- 
librium, v and p are their perturbations near equilibrium, respectively. Two dimensionless 
numbers, fluid Reynolds Re and magnetic Reynolds Rm are defined as 

n RqVq RqVq 
Re = , Rm = . 

V 1] 

In defining the normalized quantities, we use the peak driving velocity Vo, radius of the 
sphere R (a unit of length throughout the paper), plasma mass density po, the kinematic 
viscosity v and the magnetic diffusivity r\ (all three assumed to be constant and uniform). 
Eq. (1) is the Navier-Stokes equation describing the equilibrium velocity v. Since we are 
interested only in a linear (kinematic) stage of dynamo, we do not include the Lorentz force 
due to the dynamo field in Eq. (1). Eq. (2) is the Navier-Stokes equation linearized near the 
equilibrium velocity v. It constitutes an eigenvalue problem for the velocity perturbation 
v and its growth rate j v . By solving Eq. (2), one can establish the hydrodynamic stability 
properties of the equilibrium velocity v. Eq. (3) is the kinematic dynamo problem for the 
unknown dynamo magnetic field B and its growth rate 7&. Note that the growth rates in 
Eqs. (2) and (3) are normalized by the corresponding inverse diffusion times: 7^ is given in 
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units of v jR\ and 7b is given in units of r\jR\. Normalization of magnetic field B is arbitrary- 
due to linearity of Eq. (3) with respect to B. 

We restrict our study to the axisymmetric equilibrium flows only, which have dependences 
of the form v(r, 0) in the spherical system of coordinates (r, 9, <fi). Here r is the normalized 
radius (0 < r < 1), 9 is the polar angle (0 < 9 < it) and <fi is the azimuthal or toroidal 
angle (0 < < 2ir). Exploiting the geometry of the problem, we expand the divergence-free 
fields v, v and B in a spherical harmonic basis [26] and substitute these expansions into 
Eqs. (l)-(3). The resulting equations and methods of their numerical solution are given in 
Appendix A. 

Eqs. (l)-(3) are supplemented by appropriate boundary conditions. In our model, the 
boundary condition for equilibrium velocity v is specified by the driving velocity profile at 
the sphere's wall: 

V L = l=^W e ^ 0<#<7T, (5) 

where v^(9) is a function of the polar angle 9 with physical restriction t^(0) = v^ijx) = 0. 
In the present study we consider only the flows of von Karman type, i.e., the flows, whose 
driving azimuthal velocity v^iO) is antisymmetric with respect to equator {9 = ir/2). In 
such flows, the dynamo growth rate 75 is purely real near the dynamo onset, which is shown 
in Sec. IV. As a result, the critical magnetic Reynolds number Rm cr can be easily found by 
setting 7 b = in Eq. (3). For numerical convenience, we use the Fourier expansion of v,p{9)\ 

v <t>(@) — a « sm n ^ = a 2 ^ sin 2^ + /3 4 sin 4^ + /3 6 sin 6^ + . . . ^ , max tv>(#) = 1- (6) 

even n 

Here we keep only even harmonics of 9 due to equatorial antisymmetry of function v<f,(9), 
and introduce parameters f3 n = a n /ci2- In addition, by adjusting coefficient 02 we normalize 
t>^(#) so that its maximum is 1. With this definition of driving velocity t></,(#), the result- 
ing equilibrium flow (as well as corresponding Rm cr ) is uniquely determined by a set of 
independent parameters (3 n and fluid Reynolds number Re. 

The boundary condition for velocity perturbation v in case of impenetrable, no-slip wall 

is 

vUx = 0. (7) 

Specification of boundary conditions for magnetic field B plays an essential role in dynamo 
studies. As shown in Ref. [33] for the model relevant to MPDX, the critical magnetic 
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Reynolds number Rm cr is sensitive to the wall magnetic permeability but not affected by 
the wall resistivity In our present study, we assume the non-ferritic insulating wall, so 
that the normal component of the electric current is zero at the boundary, and the normal 
component of the magnetic field matches the external vacuum solution. These conditions, 
along with Eqs. (5) and (7), can be conveniently represented in terms of a spherical harmonic 
basis. The corresponding equations are given in Appendix A. 

We briefly describe the optimization procedure used in the study. First, we solve Eq. (1) 
to find the axisymmetric equilibrium velocity field v for a given set of driving parameters (3 n 
and fluid Reynolds number Re. Then, we solve Eq. (2) using v to find the eigenvalue with 
the largest real part Re{7„}. Depending on this eigenvalue, the equilibrium flow is classified 
as stable (Re{7„} < 0), marginally stable (Re{7t,} = 0), or unstable (Re{7„} > 0). The 
hydrodynamically unstable flows are not examined in our kinematic dynamo study, since 
they will develop into non-axisymmetric structures, which contradicts our analysis. We note 
that the instability of these flows is usually due to the modes with azimuthal numbers m = 1 
or m = 2. In addition, the growth rate j v is always real near the instability threshold. 

As the third step, we consider kinematic dynamo problem given by Eq. (3) with velocity 
v. Our calculations show that the fastest growing (or least decaying) dynamo modes have 
azimuthal mode number m = 1 (generally corresponding to equatorial dipoles), so we restrict 
our consideration to these modes only. In addition, for the equilibrium flows of von Karman 
type (with equatorially antisymmetric azimuthal velocity) the dynamo growth rate % is 
always real near the dynamo instability threshold, so when Rm = Rm cr the growth rate is 
zero 7f, = 0. This circumstance allows us to significantly simplify the procedure of finding 
Rm cr . We solve Eq. (3) with 7^ = as a generalized eigenvalue problem for Rm cr . The 
minimal positive eigenvalue found corresponds to the required Rm cr . If there are no positive 
numbers among calculated eigenvalues, then the dynamo cannot be excited or it is excited 
with Im{7 fe } 7^ 0. 

These three steps define Rm cr implicitly as a function of fluid Reynolds number Re and 
independent driving parameters ft n . For a given Re we search for a minimum of Rm cr in 
a multi-dimensional space of the driving parameters (5 n and determine the optimized flow. 
Multi-dimensional numerical minimization in our study is realized via the downhill simplex 
method (also known as Nelder-Mead or amoeba method [34]). 
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FIG. 2: Profiles of driving boundary velocity iu(#) for flows I, II and III. 




FIG. 3: Equilibrium structures of flows I, II and III for fluid Reynolds number Re = 150. The left 
half of each figure shows stream lines of poloidal velocity v po i superimposed on its absolute values 
depicted in colors, the right half shows a contour plot of azimuthal velocity (dashed curves 
denote values of iu < 0). The vertical central lines represent the axis of symmetry. 

III. DYNAMO ONSET 

In this section we study the effects that influence the dynamo onset. For this purpose we 
consider three axisymmetric equilibrium flows (denoted as flow I, II and III, respectively) 
driven by different profiles of boundary azimuthal velocity v$(9) as shown in Fig. 2. The 
corresponding equilibrium flow structures for Re = 150 are presented in Fig. 3. Flow I is 
driven by boundary velocity v^O) = sin 2$, i.e., only the first harmonic in Eq. (6) is kept. 
Driving velocity profile of flow II consists of half-sines with opposite signs in the intervals 
< 9 < 7i / 3 and 2tt/3 < 9 < tt. Boundary velocity of flow III is taken from Ref. [10]; it is 
similar to that of flow II, but has several reversals (changes of sign) in equatorial region. The 
driving Fourier coefficients for these flows are given in Table II. The coefficients are truncated 
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TABLE II: Fourier coefficients a n of driving velocities iu(0) for flows I, II and III. 
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FIG. 4: Critical magnetic Reynolds Rm cr as a function of fluid Reynolds Re for flow III. 

at n — 20, this is the number of harmonics used in the kinematic dynamo calculations. We 
note that among these flows only flow III results in dynamo action (dependence of critical 
magnetic Reynolds number Rm cr on fluid Reynolds number Re is given in Fig. 4). The 
reasons for this are analyzed below. 

As established in Refs. [25, 26], the pure toroidal flow cannot excite a dynamo in a 
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FIG. 5: Ratio of poloidal to toroidal kinetic energies E po i/Et or as a function of fluid Reynolds 
number Re for flows I, II and III. 

spherical geometry, a finite poloidal component of flow is necessary too. Therefore, dynamo 
onset should depend on the ratio of poloidal to toroidal flow amplitudes or, equivalently, 
the ratio of poloidal to toroidal kinetic energies. For the counter-rotating axisymmetric flow 
of Dudley and James [4] it was found that this ratio must be in a certain range in order 
to excite a dynamo. In a boundary driven flow, this ratio is determined by the details of 
the driving velocity profile and value of fluid Reynolds number Re. Fig. 5 shows the ratio 
of poloidal to toroidal kinetic energies E po i/E t0 r as a function of Re for flows I, II and III. 
The figure suggests that the failure of flow I to support a dynamo most likely results from 
insufficient relative intensity of its poloidal component. 

Fig. 5 is obtained assuming that flows I, II and III are axisymmetric at 100 < Re < 300. 
However, in reality their axial symmetry breaks when they become hydrodynamically un- 
stable. The corresponding threshold values of fluid Reynolds number, above which the non- 
axisymmetric instabilities develop in these flows, are Rej ~ 115 (instability with azimuthal 
mode number m = 2), Ren ~ 207 (m = 2) and Rem ~ 305 (m = 1). 

The presence of poloidal flow component with large enough amplitude compared to 
toroidal component is a necessary but not sufficient condition for dynamo onset. The simi- 
lar driving velocity profiles of flows II and III produce similar ratios of poloidal to toroidal 
kinetic energies E po i/E t0 r (Fig- 5). However, they are completely different from the dynamo 
point of view: flow III leads to a dynamo action, while flow II does not. Such difference can 
be explained by the details of driving velocity profiles. Namely, the presence of reversals in 
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v<j,(9) in equatorial region of flow III appears to be crucial for the dynamo onset. All the 
optimized flows found in Sec. IV possess this property. Currently, the role of these reversals 
for dynamo onset is not understood completely. 



IV. RESULTS OF FLOW OPTIMIZATION 

In this section we report the results of the flow optimization, which is performed to 
minimize the critical magnetic Reynolds number Rm cr required for the dynamo onset. First, 
we consider equilibrium flows driven by f </,(#) with two lowest Fourier harmonics in Eq. (6): 



Such equilibrium flows are uniquely determined by two independent parameters: Fourier 
coefficients ratio /3 4 = 04/02 and fluid Reynolds number Re. The corresponding Rm cr is 
also a function of these two parameters. The contour plot of this function Rm cr (Re, j3 4 ) is 
shown in Fig. 6. Dynamo action exists only in a bounded domain of the parameter space, 
approximately in the range 100 < Re < 270 and 0.7 < f3 4 < 1. Note that if f3 4 = 0, then 
V 4>(Q) ~^ sin 2$, and the dynamo action is not possible (flow I from Sec. III). Increase in fluid 
Reynolds number Re makes the background equilibrium flow hydrodynamically unstable 
with respect to non-axisymmetric modes (shaded area), and our kinematic dynamo analysis 
is not valid in this case. Scanning Re we determine values of (3 4 , which lead to optimized 
equilibrium flows minimizing Rm cr ; these values form a solid black curve in Fig. 6. For 
Re > 245 the optimized flows are at the boundary of hydrodynamic stability (dashed curve). 
The global minimum of the critical magnetic Reynolds number in this case is Rm cr pa 237 
achieved at Re ~ 250. 

In Fig. 6, the lower (/3 4 pa 0.7) and the upper (f3 4 pa 1) dynamo/no dynamo transitions 
have different behavior of Rm cr . This difference is clarified in Fig. 7, which shows curves of 
7 6 (i?m) at Re = 150 for several values of f3 4 . Rm cr increases as f3 4 decreases from 0.75 to 
0.7. This corresponds roughly to tilting the curve ^(Rm) to the right in Fig. 7. When f3 4 
approaches 0.7 from above, Rm cr goes to infinity, and there is no dynamo for f3 4 < 0.7 (the 
lower dynamo/no dynamo transition). Increasing f3 4 from 0.9 to 1.0 corresponds roughly to 
shifting the parabolic curve ^(Rm) down, below the 75 = line in Fig. 7. This explains 
why the upper dynamo/no dynamo transition at j3 4 ~ 1 occurs at finite values of Rm cr . 




max v<f,{6) = 1. 
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FIG. 6: Contour plot of critical magnetic Reynolds number Rm cr as a function of fluid Reynolds 
number Re and driving parameter (3^ = 04/02- Contours of Rm cr are shown. The shaded area de- 
notes the hydrodynamically unstable region, with stability boundaries shown for azimuthal modes 
m = 1,2,3 (curves labeled with symbols). Each point of the solid black curve corresponds to the 
optimized stable flow that minimizes Rm cr at a given value of Re. Points of the dashed black 
curve correspond to the optimized flows at the boundary of hydrodynamic stability. Symbol "x" 
denotes the point (Rm cr ~ 237, Re ~ 250), at which the global minimum of Rm cr is achieved. The 
segmentation of the dynamo/no dynamo boundary is due to discrete scan of the plane Re — (3^. 

We perform similar analysis for the flows with different number of Fourier harmonics in 
driving boundary velocity v^{9). The respective cases are marked according to the number 
of the highest non-zero harmonic in Eq. (6) (for example, N = 8 means the flows driven 
by all even harmonics up to sin 86). The results of the analysis are summarized in Figs. 8 
and 9 for six cases with values of N ranging from iV = 4 to N = 14. The curves in Fig. 8 
show the lowest possible Rm cr achievable for a given Re by optimizing the flows with a 
different number of driving harmonics. In all cases, Rm cr decreases with increasing Re for 
the optimized hydrodynamically stable flows (solid curves) and reaches a minimum for the 
flows at the stability boundary (dashed curves). The overall lowest Rm cr « 200 is obtained 
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FIG. 7: Real part of dynamo growth rate Re{7(,} as a function of magnetic Reynolds number Rm 
for fluid Reynolds number Re = 150 and different values of driving parameters /?4 = 0.7 (solid 
curve), /?4 = 0.75 (dashed curve), (3^ = 0.9 (dashed-dotted curve) and (3^ = 1.0 (dotted curve). 

at Re ~ 240 for N = 14. Fig. 9 shows the dependences of driving Fourier coefficients a n 
on fluid Reynolds number Re in the optimized flows. Note that only two coefficients a2 and 
CI4 are relatively large in all cases (04 ~ a§ ~ 0.6), the magnitude of others is normally less 
than 0.2. 

Fig. 10 demonstrates the samples of optimized profiles of boundary velocity v ^(8) for dif- 
ferent number of driving harmonics (optimization is done at Re = 150). The corresponding 
structures of equilibrium flows are shown in Fig. 11. The optimized driving velocities v^O) 
exhibit changes of sign in equatorial region. As mentioned in Sec. Ill, the presence of such 
reversals appears to be crucial for the dynamo action. Also we note that for the cases with 
N = 10, N = 12 and iV = 14, both the optimized profiles of v^(8) and the corresponding 
flow structures are very similar. This indicates that at this stage, the optimized flows are 
not strongly affected by higher driving harmonics. 

Fig. 12 shows the dependences of dynamo growth rate 75 on magnetic Reynolds number 
Rm obtained for the optimized flows from Fig. 11. A dynamo can be excited only in a 
relatively narrow range of Rm. In addition to Rm cr required for the dynamo onset, there 
is the second Rm cr 2, above which the dynamo is quenched. This is a typical indication of 
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FIG. 8: Critical magnetic Reynolds Rm cr as a function of fluid Reynolds Re for optimized flows with 
different number of driving harmonics. Solid curves correspond to the optimized hydrodynamically 
stable flows, dashed curves correspond to the optimized flows at the boundary of hydrodynamic 
stability. N denotes the number of the highest non-zero Fourier harmonic in the driving velocity 
v ( p(6). Symbol "x" denotes the overall lowest value of Rm cr « 200, which is achieved at Re « 240. 
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FIG. 9: Dependences of driving Fourier coefficients a n on fluid Reynolds Re in optimized flows 
with different number of driving harmonics. Values of Re are scanned with step ARe = 10. As 
in Fig. 8, solid curves correspond to the stable flows, dashed curves correspond to the marginally 
stable flows. 
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FIG. 10: Profiles of optimized azimuthal velocity vj,(0) at the boundary for different number of 
driving harmonics, fluid Reynolds number is Re = 150. 
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FIG. 11: Axisymmetric equilibrium flows corresponding to optimized driving velocities from Fig. 
10 at fluid Reynolds number Re = 150. Notations are the same as in Fig. 3. 
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FIG. 12: Dependences of real (solid curves) and imaginary (dashed curves) parts of dynamo growth 
rate 75 on magnetic Reynolds number Rm for the optimized flows shown in Fig. 11. Calculations 
are done for the fastest dynamo mode (with azimuthal mode number m = 1). 

slow dynamos [35]. The structures of the growing dynamo eigenmodes at Rm = 400 are 
presented in Fig. 13. The excited dynamo field outside the sphere has dipole-like structure 
and the axis of dipole is perpendicular to the axis of flow symmetry. 



V. SUMMARY 



We have numerically found the optimized laminar axisymmetric flows in a sphere that 
minimize the critical magnetic Reynolds number Rm cr required for the dynamo action. The 
flows are solutions to the Navier-Stokes equation with von Karman type boundary conditions 
specified by azimuthal velocity profiles v^(9) at the sphere's wall. In this class of flows, the 
overall minimum of Rm cr ~ 200 is obtained at the fluid Reynolds number Re ~ 240 when 
the flow is hydrodynamically marginally stable. Solving the kinematic dynamo problem for 
the optimized flows, we have determined that the dynamo action is quenched above the 
second critical magnetic Reynolds number Rm cr 2, which is typical for slow dynamos. 
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FIG. 13: Magnetic field lines of the fastest kinematic dynamo eigen-modes obtained at Rm = 400 
for the optimized flows shown in Fig. 11. Thickness of the lines is proportional to the magnitude 
of the field. Vertical lines represent the axis of symmetry of the flows. Horizontal lines (shown in 
darker color) denote the axis of these equatorial dipole-like dynamo fields. 

We have shown that in a boundary driven flow the dynamo can be excited only when the 
ratio of poloidal to toroidal kinetic energies in that flow is sufficiently high. This condition 
is necessary for the dynamo onset but does not guarantee it. In addition, the dynamo 
onset is very sensitive to the details of driving velocity profile. In all cases explored here, 
axisymmetric flows of von Karman type sustain dynamo only if the corresponding driving 
velocities have reversals (sign changes) near the equator. The effect of these reversals on 
dynamo onset is not fully comprehended. 

Our results suggest that the dynamo excitation can be demonstrated in an experiment 
with controllable laminar plasma flows, such as MPDX. Simple estimates show that the 
dynamo regime with fluid Reynolds number of Re = 150 and magnetic Reynolds number 
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of Re = 400 can be reached in an argon plasma with peak driving velocity of Vo = 5 km/s, 
density n = 10 18 m~ 3 , electron and ion temperatures of T e = 10 eV and Tj = 1 eV, 
respectively. These parameters will soon be achievable in the MPDX. 
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Appendix A: Equations and numerical methods 

In a spherical harmonic basis the divergence-free vector fields v, v and b are represented 



as [26] 




(Ala) 



i=i 




(Alb) 



i=i 




(Ale) 



i=i 




(A2a) 



l=m 




l=m 



(A2b) 




l=m 





(A3a) 



l=m 




l=m 



(A3b) 




L b r 



rn 



(A3c) 



l=m 
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where Y™ are the spherical harmonics defined in terms of the associated Legendre polynomi- 
als as Y l m (6, 4>) = Pl n (cos6)e im< ^. Since the equilibrium velocity field v is axisymmetric, and 
Eqs. (2), (3) are linear in v and B, each azimuthal mode m of v and B can be considered 
separately. The summations in Eqs. (Al)-(A3) are truncated at some (generally different) 
spherical harmonics of degrees Lq, L v and Lb, respectively. In numerical calculations we 
discretize the unknown functions of radius r (0 < r < 1) on a uniform radial grid with N r 
intervals and apply the finite difference method of the second order. All results reported 
in the paper are obtained with L = L v = L b = 20 and N r = 50. This spatial resolution 
appears to be sufficient for the cases under consideration. The convergence of the numerical 
schemes is checked by comparing simulations at different resolutions. 

In the following subsections we consider in more details Eqs. (l)-(3) for the fields given 
by Eqs. (Al)-(A3) and methods of their solution. 



1. Equilibrium velocity 



Substituting Eqs. (Al) into Eq. (1) and using the orthogonality properties of spherical 
harmonics, we obtain for 1 < I < L : 



A? a , = ReA^Y, 

3=1 k=l 



d ( SjA k s k + tjt k \ ( A k s k dsj t k dtj 



dr \ r 2 



jik 



3=1 k=l 

Here, A ; is the differential operator 



dtk _ 9sfc 
3 dr 1 dr 



(A4a) 
(A4b) 



x= d 2 1(1 + 1) 
dr 2 r 2 ' 

A™ is the numerical factor (we give the general expression for arbitrary m, since it is used 
below) 



AT 



(2/ + l)(/-m)! 
21(1 + T){1 + m)V 



and elements C\ 3k are defined as 



/BY f)Y° 
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The boundary conditions for functions S/(r) and ti(r) in Eqs. (A4) follow from the absence 
of a singularity in the velocity field at the center of the sphere and Eq. (5): 

= 0, *i| r=0 = 0, (A5a) 

r=0 

= 0, Ul^n. (A5b) 



si 



Sl 



r=0 Q r 

dsi 

r=l " 



r=l 

Here Tj are coefficients of the expansion of driving velocity v<j,{6) in terms of (—dY [ °/d9): 

L ° flyO 

= "I>-gp ( A6 ) 

Coefficients 77 are uniquely determined by Fourier coefficients a n from Eq. (6). Indeed, 
comparing Eqs. (A6) and (6) we have 

1=1 n=l 

where for generality we assumed Fourier expansion of v$(9), which includes both even and 
odd harmonics of 9. Multiplying both sides of this equality by (— A® sin6> dY^/dO), integrat- 
ing over < 6 < n and using orthogonality properties of the spherical harmonics we arrive 
at invertible matrix transform: 

n = ^ F ln a n , 1 < I, n < L 

n=l 

with matrix elements F\ n given by 

21 + 1 f dY° 
Fin = —7^r, — 7T / sin n9 sin 9—^-d9. 



21(1 + 1) J 09 
o 

In order to solve Eqs. (A4) with boundary conditions given by Eqs. (A5), we use an 
iterative scheme. The iterations are organized in the following way. First, by inverting the 
operators and A/ (which are tridiagonal square matrices in finite difference representa- 
tion) we bring Eqs. (A4) to the form 

x = f(x), x=(sj,ij), 

where f(x) denotes the nonlinear right-hand side of Eqs. (A4). Then we construct the 
iteration step: 

x (p+i) = x (p) + a [ f ( x (p)) _ x (p)] ; ( A7 ) 
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where a is a constant, chosen to guarantee convergence of iterations. The iterations can be 
initialized with some profiles of s\°\r) and t\°\r) satisfying the boundary conditions. The 
iterations stop when 



f(x 



(p)^ 



< ex 



(p)l 



where e is the error tolerance and the norm ||x|| is defined as a sum of squares of absolute 
values of all elements in x. The results of the paper are obtained with s\°\r) = 0, t\°\r) = 
T t r l+1 i a = 0.01, e = 10~ 12 . This choice provides fast convergence to equilibrium state with 
typical number of required iterations ~ 10 3 . 

2. Hydrodynamic stability 



Substituting Eqs. (Al) and (A2) into Eq. (2), we obtain for m < I < L v : 



L v 



iPAiSj + J, 



dsj d 



jl T- -jl ^ T- ^ (iff I 



(4) US] 



AjSj + J, 



f%~<^ (A8a) 



7 (3)^_ 7 (4)~ , j(2) t ~. 



r(4)^7 , t(1)7 



Ly 



{2)Uij (l)~ ( 3 )~ .(3) 



lj Q r lj 



J ii J + J lj 3 ^ 1 lj LA 3 b ] J l 



(2)d5j , 



lj Q r ~ lj 



.(A8b) 



Here the bar above a symbol denotes its complex conjugate, 1^ ^ are functions of r deter- 
mined by equilibrium profiles Sfe(r), tk(r), 



ljj \ s k, t k ] 



i(i + i) 



Ly 

E 

fe=i 



4? } = E 



k(k + l)sfe 



fc=i 



Z(Z + l)tf fcU - M, 



r(3)r_ / i _ V k(k + l)s k 



k=l 



Ilj\ s k,tk] — 



j'O' + i) 



Ly 

E 

fe=i 



<9r 



t k M, 



kl;j 



(A9a) 
(A9b) 
(A9c) 
(A9d) 



and J^ 1 ^ are obtained from Eqs. (A9) by replacement s k — > t k and t k — > —A k Sk, i.e., 
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The elements K k ij, L k ij and M k ij in Eqs. (A9) are integrals of triple products of spherical 
harmonics: 



K, 



kij 



= J Y k Y™Yf sin OdB, 



L kij = im 



dY 



k \/-m 



09 



o 



Mi 



kij 



f)Y° (W m 



The indices in the above formulas span the following intervals: 1 < k < L , m < l,j < L v . 
Note that M k ij can be reduced to K k [j by integration by parts, namely, 

1 



Mi 



kij 



[k(k + l) + l(l + l)-j(j + l)]K klj . 



The boundary conditions for functions S/(r) and t/(r) are 

si 



\r=0 



Si 



r=l 



dsi 
dr 
dh 
dr 



r=0 



0, tA =0, 
' l \r=0 ' 



= o, f,| , = 0. 

' l \r=l 



(AlOa) 
(AlOb) 



r=l 



Eqs. (A8) and (A10) constitute an eigenvalue problem for the growth rate j v and un- 
known functions §i(r) and ti(r) describing the velocity perturbations for the azimuthal mode 
m. Applying the finite difference method, we reduce the problem to a generalized matrix 
eigenvalue equation with matrix size [2(L V — m + l)(N r — l)] 2 . This equation is solved in 
MATLAB. 

3. Kinematic dynamo 



Substituting Eqs. (Al) and (A3) into Eq. (3) we obtain for m < I < L b : 

L b 



lb S t = AM + RmA?^ 

j=m 
Lb r 



IPS 



r(2) dS 3 



lj 3 A lj Q f . 



— I — T^T- 



'j=m 



<L( T&T -i- T^ 9Sj I T^?\ 1 T W dS 3 
jl "-3 "t- q t \ A ij J-j -I- hj q t -r A lj D i J A jl q t 



(Alia) 
,(Allb) 
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where integrals 1^ ^ are defined by Eqs. (A9) with indices l,j spanning the interval m < 
l,j < L b . The boundary conditions for the non-ferritic insulating wall are obtained by 
matching the dynamo field onto a vacuum potential field solution at r = 1: 

Si\ r=0 = 0, T,| r=0 = 0, (A12a) 
-Q^ + lSij =0, T,| r=1 = 0. (Al2b) 

Formulas analogous to Eqs. (All), (A12) are originally derived in Ref. [26]. 

We transform Eqs. (All) and (A12) to a matrix eigenvalue problem by applying finite 
difference method. In this method, the boundary condition for Si in Eqs. (Al2b) is taken 
into account by using an extra (ghost) grid point to approximate the derivative at r = 1. 
The resulting matrix equation of size [(L b — m + l)(2iV r — l)] 2 is solved in MATLAB. 
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